/*******************************************************************\
| MathExt.h
|
| CREATED:		8.5.2009
| DESCRIPTION:
| 
| Extension functions to the base math library (contains more 
| complex math functions)
\********************************************************************/
#ifndef MATHEXT_H
#define MATHEXT_H

namespace Math
{

// radical inverse
inline float RadicalInverse(int p_N, int P_Base)
{
	float val = 0;
	float invBase = 1.0f / P_Base, invBi = invBase;
	while(p_N > 0) 
	{
		// Compute next digit of radical inverse
		int d_i = (p_N % P_Base);
		val += d_i * invBi;
		p_N /= P_Base;
		invBi *= invBase;
	}
	return val;
}

// returns a point on the uniform sample sphere
inline Vector3 UniformSampleSphere(float p_U1, float p_U2)
{
	using namespace Math;
	float z = 1.0f - 2.0f * p_U1;
	float r = Sqrt(Max(0.0f, 1.0f - z * z));
	float phi = 2.0f * PI * p_U2;
	float x = r * Cos(phi);
	float y = r * Sin(phi);
	return Vector3(x, y, z);
}

inline const Vector3 GetTangent(const Vector3& p_Normal)
{
	float min = INFINITY;
	int min_index = -1;

	// find the minor axis of the ray
	for(int i = 0; i < 3; ++i)
	{
		if(Abs(p_Normal.m_Xyz[i]) < min)
		{
			min = Abs(p_Normal.m_Xyz[i]);
			min_index = i;
		}
	}

	Vector3 result = Vector3::ZERO;
	if(min_index == 0)
	{
		result = Vector3(0.0f ,-p_Normal.m_Xyz[2],p_Normal.m_Xyz[1]);
	}
	else if(min_index == 1)
	{
		result = Vector3(-p_Normal.m_Xyz[2], 0.0f ,p_Normal.m_Xyz[0]);
	}
	else if(min_index == 2)
	{
		result = Vector3(-p_Normal.m_Xyz[1], p_Normal.m_Xyz[0], 0.0f);
	}

	result.Normalize();
	return result;
}

// gaussian filter function for radiance estimation
inline float GaussianFilter(const float p_SqDist, const float p_SqRadius)
{
	const float gAlpha = 1.818f;
	const float gBeta = 1.953f;

	const float gNom = 1.0f - Exp(-gBeta * p_SqDist / (2 * p_SqRadius));
	const float gDen = 1.0f - Exp(-gBeta);
	return gAlpha * (1.0f - (gNom) / gDen);
}

// cone filter function for radiance estimation
inline float ConicFilter(const float p_Dist, const float p_Radius)
{
	const float coneFilter = 1.1f;
	return 1.0f - (p_Dist / (coneFilter * p_Radius));
}

}

#endif // MATHEXT_H
